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ABSTRACT 

CL 

[T"1 , The recent discovery of a transiting short-period planet on a slightly non-circular 

orbit with a massive highly eccentric companion orbiting the star HAT-P-13 offers 
O . the possibility of probing the structure of the short-period planet. The ability to do 

this relies on the system being in a quasi-equilibrium state in the sense that the 
eccentricities are constant on the usual secular timescale (typically, a few thousand 
years), and decay on a timescale which is much longer than the age of the system. 
Since the equilibrium eccentricity is effectively a function only of observable system 
parameters and the unknown Love number of the short-period planet, the latter can 
^vq , be determined with accurate measurements of the planet's eccentricity and radius. 

However, this analysis relies on the assumption that the system is coplanar, a 
situation which seems unlikely given the high eccentricity of the outer planet. Here 
we generalize our recent analysis of this fixed-point phenomenon to mutually inclined 
systems in which the outer body dominates the total angular momentum, and show 
that (1) the fixed point of coplanar systems is replaced by a limit cycle in e?, — rj 
■ space, where is the eccentricity of the inner planet and 77 is the angle between the 

, periapse lines, with the average value of e(,, e^ av \ decreasing and its amplitude of 

variation increasing with increasing mutual inclination. This behaviour significantly 
reduces the ability to unambiguously determine the Love number of the short-period 
planet if the mutual inclination is higher than around 10°. (2) We show that for Q- 
values less than 10 6 , the HAT-P-13 system cannot have a mutual inclination between 
54 and 126 degrees because Kozai oscillations coupled with tidal dissipation would 
act to quickly move the inclination outside this range, and (3) that the behaviour of 
retrograde systems is the mirror image of that for prograde systems in the sense that 
(almost) identical limit cycles exist for a given mutual inclination and n minus this 
value. (4) We derive a relationship between e^ v \ the equilibrium radius of the short- 
period planet, its Q-value and its core mass, and show that given current estimates 
of e& and the planet radius, as well as the lower bound placed on the Q- value by the 
decay rate of e£ , the HAT-P-13 system is likely to be close to prograde coplanar, 
or have a mutual inclination between 130° and 135°. Lower rather than higher core 
masses are favoured. (5) An expression for the timescale for decay of the mutual 
inclination is derived, revealing that it evolves towards a non-zero value as long as 
&h > on a timescale which is much longer than the age of the system. (6) We 
conclude with a scattering scenario for the origin of the HAT-P-13 system and show 
that almost identical initial conditions can result in significantly different outer planet 
eccentricities, stellar obliquities and planet radii. The implications for systems with 
high stellar obliquities such as HAT-P-7 and WASP-17 are briefly discussed. 
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Table 1. Parameters of the HAT- P- 13 system 



HAT-P-13 

m, (M ) 1-22±8:S§ 
R* (Rq) 1.56 ±0.08 

5 n+ 2 -5 

°- U -0.8 

fc, 0.03 
r 9 */R* 0.076 
7? fl 2.04 
7^ de 0.04 
7* pi " [iVn = Pq] 0.06 



age (Gyr) 



+0.029 
-0.046 
+0.0006 



HAT-P-13b 

m b (Mj) 
a b (AU) 

P b (days) 2.916260 ± 0.000010 

e b 0.021 ± 0.009 

u> b 181 ± 46° 

R b (Rj) 1.280 ±0.079 

i los 83.4 ± 0.6° 

k b 0.3 

r gb /R b 0.26 

7^ de 4.83 

r circ 40(Q 6 /10 5 ) Myr 

Ta 79(Qi,/10 5 ) Gyr 

r e 15(Q 6 /10 5 ) Gyr 

r c 9200 (Q o /10 5 ) Gyr 

80(Q b /10 5 ) Gyr 

HAT-P-13c 

m c sin i ios = m™ n (Mj) 15.2 ± 1.0 

a c (AU) 1.186±g;°i| 

P c (days) 428.5 ± 3.0 

e c 0.691 ± 0.018 

ujc 176.7 ± 0.5° 



INTRODUCTION 



Transiting systems offer the opportunity to determine a wide variety of system parameters, including the mass of the transiting 
planet, its orbital eccentricity, the inclination of its orbit to the line of sight, its radius and hence mean density, and the sky 
projection o f the angle between its orbit normal and the stellar spin axis. Many other system parameters are potentially 
measurable (IWinnll2009i ). one of them being the Love number of the tra nsiting planet if th e system parameters are favourable 
( Wu fc Goldreichl 20021 ) . Th e recent discovery of the HAT-P-13 system jBakos et al.ll200^ ) provides us with such a system as 
was recently pointed out by Batygin. Bodenheimer &: Laughlin ( 20091 ). A reliable estimate of a planet's Love number in turn 
allows one to say something about the presence or otherwise of a planetary core, and hence about the formation mode; in the 
former case t he planet is likely t o have been formed via core accretion of solid material and subsequent accretion of a massive 
atmospher e JPollack et ah 1996J ). while the absence of a core supports the gravitational collapse hypothesis of giant planet 
formation l|Bosslll997l ). 

The HAT-P-13 system consists of a 0.85 Jupiter-mass planet (planet b) in an almost circular 2.9 day orbit about a 1.2 
solar-mass star, and a companion with a minimum mass of 1 5.2Mj in a 428 day highly eccentric orbit (planet c). Table [T] lists 
the relevant parameters of the system, with data taken from lBakos et al. (2009) who performed a simultaneous fit of HATNet 
and KeplerCam photometric data and Keck spectroscopic data, a process they refer to as "global" modelling. Here P b and 
P c are orbital periods, ii oa is the inclination of the orbit normal of planet c to the line of sight and ui b and uj c are periapse 
arguments. Note in particular the non-zero estimate of 0.021 ± 0.009 for the inner planet's eccentricity. The main aim of this 
paper is to explore the extent to which information can be gleaned from such a non-zero measurement, whether or not we 
know the mutual inclination of the orbits of the two planets. 

The ability to determine the Love number of HAT-P-13b relies on the system being in a quasi-equilibrium state in the 
sense that after an initial transient phase of rapid tidal evolution, the eccentricities change on a timescale much longer than 
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the age of the system and the apsidal lines of the inner and outer planets are aligned or anti-aligned. Moreover, the mass and 
eccentricity of HAT-P-13c ensure that t he quasi-equilibrium ecc entricity of the inner planet is significant and measurable. 
S ome of the theory for this is developed in IWu fc Goldreichl (|2002h for the HD83443 system, while a general theory is presented 
in Mardlinel ( 20071 ). Both of these studies, however, are for the coplanar case only. Using the HAT-P-13 system as illustration, 



we generalize the theory to non-coplanar systems in which the outer body dominates the total angular momentum and show 
that the quasi-relaxed state no longer corresponds to a fixed point in et — r\ space, where n — zu b — ro c with za b and vjc 
the longitudes of periastron, but rather corresponds to a limit cycle in this space (see, for example, Jordan & Smith 1999), 
with the average value of e b decreasing and its amplitude of variation increasing with increasing mutual inclination. As in the 
coplanar case, the relaxation timescale is around three times the circularization timescale. However, unlike the coplanar case, 
the rate of change of the argument of periastron of the inner planet plays an important role, with the limit cycle frequency 
equal to twice this quantity. In fact it is the appearance of this additional frequency which prevents the system from evolving 
to a fixed point, with terms which depend on it effectively acting as external forcing terms. Note that our analysis takes into 
account the fact that the actual mass of the outer planet in creases with incre asing mutual inclination, given the observed 
minimum mass determined via radial velocity measurements ( Bakos et al. 2009h . 

The plan for this paper is as follows. Section 2 reviews the theory for the long-term tidal evolution of coplanar systems, 
Section 3 generalizes this to mutually inclined systems, including limit cycle behaviour of prograde systems (Section 3.1), 
the Kozai regime (Section [32]), limit-cycle behaviour of retrograde systems (Section 13. 3p . the relationship between e^ v \ the 
equilibrium radius of planet b, its Q-value and its core mass fSection l3.4[) . and the timescale for decay of the mutual inclination 
in a relaxed system ( Section 13. 5[) . Section 4 presents scattering scenarios for the origin of HAT-P-13-like systems, including 
a discussion of stellar obliquity in two-planet systems CSection l4.1[) . Section 5 presents a summary, and Appendix A presents 
the orbit-averaged equations of motion for a two-planet Newtonian point-mass system up to octopole order, correct to leading 
order in the inner planet's eccentricity, the ratio of semimajor axes, and the sine of the inclination of the outer orbit relative 
to its initial orbit. 



2 LONG-TERM TIDAL EVOLUTION OF COPLANAR SYSTEMS 



In lMardlina (|2007f ) . the long-term tidal evolution of short-period planets with single companions is studied. There it is shown 
that such systems evolve on three distinct timescales, an illustration of which is given in Figures 3 and 4 of that paper. As long 
as the eccentricity of the outer planet is non-zero, the eccentricities of both planets will initially execute anti-phased secular 
oscillations until a non-zero quasi-equilib rium value is rea ched, with maxima and minima of the inner planet's eccentricity 
given by expressions (20), (27) or (28) in iMardlind (|2007l ). the choice of which depends on the system parameters, and the 
corresponding values of the outer planet's eccentricity given by expression (29 ). The equilibri um value of the eccentricity 
depends on all contributions to the rate of apsidal motion of the inner planet. In lMardling ([2007j) only the contributions from 
the ou ter planet and the post-Newtonian terms in the star's potential were taken into account, however, as lRagozzine fc Woli 
(2009) joint out, the contribution of the tidal bulge of a short-period planet like HAT-P-13b is significant and is included 



here (as are the contributions from the spin bulge of the planet and the tidal and spin bulges of the star). The equilibrium 
eccentricity for a coplanar system is given by 



>9) 



(5/4)(ob/a c 



1 - yja b /a c (m b /m c )e 



(1) 



where et,, a b , and m b are respectively the innermost planet's eccentricity, semimajor axis and mass, and e c , a c and m c are the 
corresponding values for the outer planet. Here e c = y/1 — ei and 7 
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to first-order in the inner eccentricity. Here n b is the mean motion of the inner planet, c the speed of light, R b and k b the 
radius and Love number of the inner planet, m», R*, fc* and f2* the star's mass, radius, Love number and spin frequency 
respectively, and synchronous rotation of the planet is assumed in the expression for ^ pln . The 7's are the ratios of the various 
contributions to the apsidal motion of the inner planet to the coplanar contribution of the outer planet. These are listed in 
Table[T]for the HAT-P-13 system (for m c = mj™ ln ), together with those for the star assuming a spin period of 25 days. Th e 
Love number for the star (equal to twice its apsidal motion constant) is taken to be that for an n = 3 polytrope l Sternell941 ). 
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If the angle between the apsidal lines of the two planetary orbits, 77, circulates rather than librates, the average eccentricity 
of the inner planet during the initial oscillatory phase will decrease on its own circularization timescale until its minimum 
eccentricity is zero, with the amplitude of oscillation remaining constant. This timescale is given to second order in the 
eccentricity by 



21./;,, V A-,, I \m.J\Rb) ' ^ 
where Qb is the Q-value of the inner planet. Once the minimum of the inner eccentricity reaches zero, r) will librate about 
zero or -k (the choice of which depends on the system parameters), and the oscillation amplitude will reduce to zero until the 
inner eccentricity reaches a non-zero quasi-equilibrium value. The librating phase occurs on a timescale of 2r c ir C - 

Once the oscillatory phase is over, the system will evolve to the doubly -circular state on the approximate timescale 

'-(i)(s)(sr-'«>-- <«> 

where e* is the value of e c at the beginning of this phase, and 



F(e c ) = e s c (l - y/ a b /a c (m b /m c )e c 1 + ye 3 c ) 2 = e c Al. (7) 

Note, however, that if the inner semimajor axis evolves appreciably on this timescale, ((6]) represents an upper limit only. 
Note also that r c is independent of the circularization timescale of the outer planet; if the latter is comparable to or shorter 
than t c then ((6} again represents an upper bound. For the HAT-P-13 system, we have T c j rc = 4 x 10 7 (Qj,/10 5 ) yr and 

Tc 



t c then (|6} again represents an upper bound. For the HAT-P-13 system, we have T c j rc = 4 x 10 (Q/,/10 
2.5 x 10 5 T C i r c = 10,000(Q 6 /10 5 ) Gyr, while the orbital decay timescale for planet b l|Yoder fc Peale! Il98lh 



e b 2 T C irc — 2500 Tdrc ~ 100(Q{,/10 5 ) Gyr < r c , the latter two being muc h grea ter than the age of the system. A direct 
coplanar integration using the averaged code presented in iMardling fe~Lml (|2002h (using a constant value of the radius of 



planet b) gives r a = 79 (Q/,/10 5 ) Gyr and r c = e c /e c = 9200 (Q 6 /10 5 ) Gyr, while e[ eq) decays on a timescale r e = 15 (Qb/10 5 ) 
GyrQ The latter is consistent with the estimate 

r- 1 = [1 + (87 - 4 7 ? fl ) £ g/A ]r- 1 ~ 6 r„-\ (8) 

obtained from (fj} with Aq — l + 7£c- The numerically determined timescales are listed in Table [T] together with the timescale 
for decay of the mutual inclination to its equilibrium value, n, for the case that it is initially 30° (see Section [33]). 

The estimate for r e places a lower bound on the value of Qb such that Qt/10 5 > (age/r^s)/ ln[e,,(0)/e(,(age)], where 
r ej 5 = T e (Q b — 10 5 ), eb(0) is the "initial" value of e& and ei,(age) is its value now. Given the lower bound for the age of 
the system, using r ej 5 = 15 Gyr and taking e t,(0) = 1 gives Qb > 7000, while et ( 0) = 0.1 gives Qb > 1.8 x 10 4 . Note that 



our estimate for r e is more than twice that of Batygin. Bodcnhci mer &: Lau ghlin (2009) who estimate r e ~ 6(Q;,/10 ) Gyr 



making their lower bound for Q b a factor 15/6 higher. 

In general, coplanar systems for which one may say something about the structure of the short-period planet will have 
the following characteristics: 

(1) The circularization timescale of planet b will be (considerably) less than one third the age of the system, ensuring that 
the system is sufficiently relaxed; 

(2) The timescales on which the system becomes doubly circular and the orbit decays, r c and r a respectively, will both be 
longer than the age of the system, and 

(3) The equilibrium eccentricity and the radius of the short-period planet will be measurable. 



3 INCLINED SYSTEMS 

In this Section the analysis of Mardlingl ( 2007 ) is generalized to non-zero mutual inclination. While the HAT-P-13 system 



is used for illustration, the analysis may be applied to any system for which C 1, ab/a c <g 1 and most of the angular 
momentum of the system resides in the outer orbit, the latter ensuring that the (sine of the) angle between the invariable 
plane and the outer orbit remains small. 

In order to highlight the differences between coplanar and inclined systems, we begin by presenting the results of an 
integration of a HAT-P-13-like system for which the mutual inclination is 30°, and the mass of the outer body is taken to be 
the observed minimum mass divided by the cosine of the mutual inclination. The stellar obliquity is such that the star's spin 

1 Note that IMardling fc Linl ||2002|1 code assumes a forcing-frequency-independent Q-value, a legacy of the pioneering work of 
iGoldreich fc SoteJ Jl966> whose supporting argument was based on the constancy of the Q-value over a wide range of forcing frequencies 
for the Earth, ie., for a solid body. For gaseous bodies, it seems more reasonable to use the "constant time-lag" concept originally devised 
by Darwin in which the lag angles of individual tidal components are proportional to their forcing frequencies. The latter is equivalent 
to the concept of fluid stress and its associated dissipation. Note, however, that for small eccentricities (which is the case in the present 
study), there is little difference between the formulations. 
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Figure 1. The relaxation process for an inclined system (to be compared to Figure 3 in Mardling (2007)). (a): Capture onto the limit 
cycle (red points) in ef, — r\ space. Black and green points correspond to circulation and libration respectively of rj. (b): Evolution of the 
eccentricity of planet b. (c): Detail of panel (a); the blue circle corresponds to the analytically determined limit cycle (Section [37T} . (d): 
Evolution of rj. 



axis is aligned with planet c's orbit normal. The integration is done using the "averaged" code presented in iMardling fc Lin 



(2002) which includes accelerations due to spin and tidal bulges of both the star and inner planet (both quadrupole and 
dissipative) as well as the relativistic potential of the star. Orbit-averaged expressions are used for the evolution of the inner 
orbit, while the outer orbit is integrated directly. No assumptions are made about the magnitude of any quantity except the 
ratio of semimajor axes which is assumed to be small. For now we suppress the evolution of the planetary radius and assume 
it is constant at its currently observed value. 



Figure [T] shows the analogue of Figure 3 in iMardlind 1)20071) . The initial values of r\ = rut — zu c and e^, are 70° and 0.05 
respectively. Rather than relaxing to a fixed point in et — r\ space on a timescale of 3r c i rc , the system relaxes to a limit cycle 
on the same timescale. Note that the Q-value of planet b is given an artificially low value of 50 in order to clearly demonstrate 
capture onto the limit cycle. While the modulation period is correct, the relaxation process would normally take Q&/50 as long 
as shown here. The circulatory phase is shown in black in each panel, the pre-capture libratory phase is shown in green and the 
limit-cycle phase is shown in red. Panel (c) shows an enlargement of the limit cycle centred on (e;,, r)) = (e[ a "\ 22-7r), together 
with the theoretical limit cycle derived below (blue circle). Here e^""' is the inclined-system analogue of the fixed-point value 
of et for coplanar systems, e[ e ^ (compare equations fl]) and (|19p ). While only around 1.5 relaxation timescales (ie., 4.5T c ; rc ) 
are shown here, the system was integrated for 6 relaxation timescales during which ' and it decreased by less than 0.0005 
and 0.05° respectively. Note that 20 T C i rc = 2.4(Qt / 10 5 ) Gyr, comparable to the estimated age of the system for realistic values 
of Q b . 

In order to guide the following analysis, we now present the results of several integrations of relaxed HAT-P-13-like 
systems with identical initial conditions except that the mutual inclination is varied between 0° and 50° (higher inclinations 
will be discussed in Sections 13.21 and I3.3|l . The initial value of is taken as e^ v ^ , and the apsidal lines are taken to be 
aligned. Initial angles are measured with respect to the initial orbit of planet c, that is, the relative inclination and longitude 
of planet b are specified, with the zero in longitude coinciding with the apse of planet c. Since planet c's orbit contains 98% 
of the total angular momentum, it coincid es approximately with the invariab l e plan e and there is very little change in its 
inclination as Figure |3je) shows. Following Batygin. Bodenheimer fc Laughlin ( 20091 ) we take the Love number of planet b 
to be 0.3, representative of a range of planetary structures according to their Table 1, and its radius of gyration, r g t, to be 
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Figure 2. The dependence of the relaxed state on mutual inclination. Note the different scales used for e;, in panels (a) and (b). (a): The 
solid curves represent numerically integrated solutions for % = 0°, 10°, 20° and 30°, while the red dashed lines are the instantaneous 
theoretical average values of e;,, e^ av \ given by equation (1 1 Q D . The artificially low Q- value allows the transient behaviour to die away 
quickly. Note that the modulation frequency and amplitude are independent of Qj, and are given by 2W U and A/2W U , with W M and 
A defined in equations d 1 6 D and J 1 5 1) respectively, (b): Integrated solutions for ib = 40° and 50° (black curves) together with ib = 30° 
for comparison (blue curve). The cases i b = 30° and 40° have similar averages, while the amplitude of variation of the i b = 50° system 
is relatively large. Its decay is associated with the relatively rapid decay of the semimajor axis of planet b for that case. Also shown is 
for if, = 40° (red dotted line); clearly equation I I19II is not reliable for this case. 
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Figure 3. Variation of the elements for if, = 0°, 10°, 20° and 30°. Colours in panel (b) are such that black, red, pink and blue correspond 
to 0°, 10°, 20° and 30° respectively. See text for a discussion of each panel. 



0.26-Rb, appropriate to an n = 1 polytrope. For this set of experiments, an artificially low Q-value of 10 is used to hasten the 
evolution towards the relaxed state. The stellar obliquity relative to the invariable plane normal, f?*, is set to zero initially, 
and the stellar spin period is 25 days. The stellar Love number and radius of gyration r a » are taken to be 0.03 and 0.076-R, 
respectively, the latter appropriate to an n — 3 polytrope. 

Figure [2ja) compares the evolution of e& for it, = 0°, 10°, 20° and 30°. The two main features of this plot are that 
the average eccentricity decreases with increasing mutual inclination, while the amplitude of its oscillation increases. Both 
these quantities can be determined from the analysis below, as can the modulation period of ej. Note that all three are 
independent of the Q- value of the planet (for small sin 77 and constant planet radius) so that Figure [2] represents their true 
values. Figure [2jb) shows the evolution of eb for ib = 40° and 50° (black curves) as well as for it, = 30° for comparison (blue 
curve) . 

Figure[3]shows the evolution of the other orbital elements for each value of ib - Panel (a) shows the evolution of the argument 
of periastron, each offset by 360° for clarity. Note in particular that for these initial inclinations, u>b is approximately constant, 
and is given by (|16[1 . Panel (b) shows the libratory behaviour of r/. Unlike the coplanar case which evolves to a constant value 
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such that sin 77 = —WT/\W q \ (and hence is proportional to Q^ 1 ), the quasi-relaxed state for systems with non-zero mutual 
inclination is oscillatory (although the offset is still Sin _1 [— Wr/|W^|]). The amplitude of oscillation is given by A/2W U (times 
180/71"), with A and W u defined in equations (|15|) and (|16[) respectively. 

Panel (c) shows the evolution of the longitude of the ascending node of planet b for all inclinations. Note that while planet 
b's precession rate is proportional to m c cosib (equation l)A15|l ). m c = m™ !n /| cos ib\ with m™ n held constant so that the 
precession rates are all identical. The precession period is 7281 years. Panels (d) and (e) show the evolution of the inclinations 
of orbits b and c relative to the invariable plane, ib, and i c respectively, justifying our assumption in the analysis below that 
they are constant and that sini c <C 1 (see Section [4.11 for analysis of the influence on ib of non-zero #*). Not shown is the 
variation in planet c's eccentricity; this is effectively constant with an amplitude of variation of 0.003. 

Finally, panel (f) shows the stellar obliquity relative to orbit b (the angle between the star's spin axis and planet b's orbit 
normal), ifr*b- In Section [4.11 we show that the maxima and minima of are \ib ±0*|. In fact, varies (the star nutates) 
due to torques between the misaligned stellar spin bulge and b's orbit, consistent with the variations seen in panel (f). 

In light of th e recent discove ries of high (sky-proj ected) stellar obliquities in the systems HAT-P-7 I Narita et al. 20091 : 



In lignt or tne recent discoveries ot nign ( sky-projected ) stellar obliquities m tne systems tiAl-f-i l|l\anta et aLIMUU 
Winn et al-lboosl l and WASP-17 (I Anderson et alJl201Ch . we discuss further the dynamics of stellar obliquity in Section |4~T 



3.1 Limit-cycle behaviour of prograde orbits 



In Appendix A we give the orbit-averaged disturbing function for a two-planet Newtonian point-mass system up to octopole 
order for arbitrary planet b inclination, and from this derive the equations governing the secular evolution of the elements in 
the absence of perturbations. Only leading order terms in et,, a,b/a c and sini c are retained, each of which are of order 0.01 
for the HAT-P-13 system whether or not the reference plane is taken as the initial plane of planet c or the invariable plane. 
Equations (|A12I) . (|A14[) . (|A16[) and (|A18[) should be compared with equations (4)-(7) of iMardlingl ( 20071 ) for the coplanar 
case, noting that here the inclination functions are such that f n (0) — 1, n — 1, 2, 3 and (72(0) = 0. 

Our aim now is to write down equations governing the dominant long-term behaviour of the two-planet system under 
the action of tidal dissipation, spin-orbit coupling and the relativistic potential of the star. In particular, we wish to study the 
behaviour of the system in the eb — r\ plane in order to determine whether a relaxed system with non-zero mutual inclination 
is able to tell us anything about the internal structure of planet b. Noting from Figure [3] (and the following analysis) that 
Cob = iz>b — Clb is approximately constant, and that as long as r\ librates, the argument £ = vjb + vo c — 2Qb = 2uib — 77 so that 
C, — 2ujb, the equations governing e& and 77 are approximately 



h = — [Wt e b + Wo e c sin 7;] + A e b sin(2dj b t) + B e c sin(2w 6 /j), 
and 



7? = 
with 



W q - Wo 



of—) cos 7/ + A cos(2ujbt) +£?( — ] cos(2tJbt) 
\ebJ J \ebJ 



(9) 
(10) 



LOb — 7) + ~&Jc — (lb — t&c 



(11) 



and we have taken the time origin to coincide with Wf, = 0. Here Wt = T c i rc is the inverse of the tidal circularization timescale, 
and both W q and W a have analogues in the coplanar theorj@ and are such that 



W q 

and 
W, 
with 

W(l : 



t I- \ f m b\ a-b -1 t i- \ I 3 
h{ib)-\ — )J—e c -hM + lSo 
\m c J V a c 



Wn = A • Wn 



o =HSK 2/2(ib)H/n ' 



(12) 



(13) 



(14) 



equal to minus the precession frequency divided by cosij, (see equation (|A15I) V Parameters which are zero for coplanar systems 



A=% sin 2 i b Wn and B = j (—) e~ 2 g 2 (ib)Wn, 

4 \ CLc / 



(15) 



while 



2 The subscripts q and o stand for quadrupole and octopole respectively. 
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U) b : 



COS If, 



(m b \ fab -i t ,. , 
\m c J V a c 



Wn 



(16) 



does not appear explicitly in the coplanar theory. Note that since ab/a c , it,, i c and e c are all approximately constant over a 
few times the circularization timescale (at least for small e b ), each of Wt, W q , W , Wq, uJb, A and B are also approximately 
constant on this timescale. 

Referring to Figure [2] showing the relaxed behaviour of e b for different ib, our aim now is to characterize the limit 
cycle demonstrated in Figure [1] by determining the average value of e b as well as its amplitude as a function of the mutual 
inclination. Equations ([9]) and (|10[) together represent a damped nonlinear non- homogeneous system of ordinary differential 
equations with pe riodic coefficients for wh ich the existence of periodic (limit-cycle) solutions are suggested by Figures [2] and [3] 
(see, for example, Jordan fc Smith ( 19991 ) for a discussion of such systems). Assuming a limit-cycle frequency of 2u> b and zero 



phased that r\ librates rather than circulates (this is true for ib <, 33° for the HAT-P-13 system as determined from numerical 
integrations), and that relaxed values for e b and r] are 7r/2 out of phase (a reasonable assumption given the form of ([9]) and 
(jTOj) ) , the amplitude of variation of e b and r\ on the limit cycle can be estimated by putting 



4 iC V) = 4^ l 1 - A ^ cos(2w(,t)] and r/ (ic) (t) = Vav + A tc sin(2w 6 t), 



(17) 



where e b and r\ av are defined in the following procedure. Substituting these expressions into ((9} and (I10[) . assuming Ai c <C 1 
and r)( lc ' <C 1, matching constant terms and those with phase 2&>t and neglecting terms with phase 4wt, we obtain 

A + Be B /4 av) A + (g 2 /f 2 )W q _ | sin 2 i„ + fa/f 2 )A 



Aic = 



2dj b + W e c /e b av) 



2dj b + W q 2 [cos i b + y/ab/acimb/mjsc- 1 fi(i b )} + A ' 



(18) 



where A is defined in (|12p , together with 



= (W /W q 



(5/4)(ob/a c ) e c e c 2 ■ f2(i b ) 



fa(ib) - \/a b /a c (m b /m c )e c 1 • fi(i b ) + 'ys't 



which reduces to ([T]) when ib = 0, and 

el 



Vav = -W T /W q = -l Jb 



5 Q b A 



(19) 



(20) 



Expressions (|19|) and (|20|l should be compared with equations (36) and (46) in lMardlinei 1 20071 ) . The minimum and maximum 



values of e b are therefore 



mtn,max 



(1 ± Al, 



(21) 



In fact the true amplitude tends to be around ^Ai c , perhaps due to the neglect of terms proportional to cos(4Wf,) and sin(4a>(,) 
(these are associated with the shape of the true limit cycle - see Figure [1]). We therefore replace Ai c with 



Al 



3 ^lc 



(22) 



in the analysis and figures that follow (as well as Figure [T]). Figure |4j a) plots e^ a "' from (|19[) (black solid curve) together with 



er™ and e^ ax from (J2TJ) (black dashed curves) for i b < 33° 



the range of values of ib corresponding to libration of 



i] around r\ = r/ av . Also plotted are maximum, minimum and average values of e& from a series of numerical integrations for 
HAT-P-13-like systems with 10° ^ i b ^ 50° (red solid a nd dot-dashed cu rves). The value i b — i b ' corresponds to A* c = 1 
(compare with the discussion following equation (23) in Mardline (2007)). For ij^ 1 ' <, it < i^ C2 ' ~ 46°, r\ circulates, these 
values of i b corresponding to |-4* c | > 1, while for i^ 2 ' < i b < i^ — 54°, rj again librates, this time around r\ — ir + r\ av . Here 
iff corresponds to the minimum value of i b for which Kozai oscillations occur for a given value of 7 (see Section f3 . 2 [I . 

While it is straight forward to dete rmine e"" n and e b nax for all possible (non-relaxed) coplanar configurations and their 
associated fixed points ((Marching 2007]), it not clear how to do this for inclined systems and their associated limit cycles for 



systems with i b > i b > that is, systems which r/ either circulates or librates around n. Here we content ourselves with trial 
and error expressions, obtained by trying different integer values of m and n 2 in e b av \n\Aic + n 2 ), and comparing these 
with numerical solutions. The results are summarized in Table [2] and plotted in Figure |4j a), the latter (blue dashed curves) 
showing good agreement with numerical solutions (red dashed curves). The angle if ~ 39° corresponds to the point where 



(av) 



and A\ c change sign and as a consequence, the e. 



r™ and e, 



max 
b 



curves cross. Given the form of these expressions and the 



accuracy with which they fit the numerical data, it seems likely that they are generic. 



Figure |4jb) demonstrates the dependence of e 



(av) 



and e™ ax on the Love number of planet b. Note that for the case 



k b — we have also set fc* to zero (although in the k b = 0.3 case the quadrupole moment of the star contributes only 1.4% to 
the total value of 7, assuming the star spins with the same period as the Sun). We conclude that the analysis and conclusions 



One can alternatively leave the frequency and phase as parameters to be determined using the procedure described here. 
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min (av) max min (av) max 

6 6 6 ' _ : ~ 




I 8 

i b (deg) i b (deg) 

Figure 4. (a): Comparison of theoretical (black curves) and numerical values (red curves) of ' (solid curves), e™ ln and e™ ax (dashed 
curves for theoretical and dot-dashed curves for numerical) as functions of if,. The mass of planet c is m™ In / cosi;,. Note that we have 
scaled the expression for Ai c given by II18I I by a factor of two thirds to improve the fit. The theoretical estimates for e^"" and e™ m are 

listed in Table[2] while that for ejj ^ is given by Jl9j> . Note that 77 librates around zero for ^ % < ij. , circulates for i[ C1 ' < it < i^ 2 ^ 
and librates about tt for it > K 2 \ (b): Theoretical estimates for the dependence of e^ av \ib), e™ in {ib) and e™ ax (if,) on kf,. Since the 
amplitude of variation of e;, in the relaxed state is small for ib < 10°, a measurement of ej gives a fairly accurate estimate of fc;, for those 
inclinations. 



Table 2. Eccentricity maxima and minima. 
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max 
e b 
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< Al < 1 
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e ^ 
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A* c > 1 
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A*c < - 1 


circulation 




-e ( r\Al + l) 


2 4^(^-2) 


„H 


< 


-1 < A* lc < 


libration around n 



o nBatvoin, Bodenheimer & Lauahlin (200 A) are valid as long as the mutual inclination of planets b and c is less than around 
10°; for higher values of ib, a measurement of eb does not unambiguously determine kb, although one can make arguments 
about the likelyhood of a system being near the top or bottom of the modulation cycle of et. While a Rossiter- McLaughlin 
estimate of the stellar obliquity will not help to constrain ib (see discussion in Section ^. it should be possible to fit combined 
transit and spectroscopic data for the mutual inclination IjNesvorn vl 12009). 

We now argue that for realistic values of Qb, the HAT-P-13 system cannot have a mutual inclination in the range 54° to 
126°. 



3.2 The Kozai regime 

Figure [5] shows the evolution of eb, ab, ib and r\ = zub — for the case ib(t — 0) = 60°, et,(0) = 0, af,(0) = 0.05 and Qb = 100 
(all other parameters are the same as for the cases presented in Section [3]). The system exhibits Kozai oscillations, detail of 
which is shown in panels (b) and (e) for et and ib respectively. However, given that there is no commensurate variation of ab, 
these come at the price of severe tidal decay of the orbit of planet b as shown in panel (c) . An upper bound for the timescale 
for decay of the semimajor axis of planet b's orbit due to planetary tides is r a — T C i rc je\, where r C i rc is given by |[5j. With 
eb at least an order of magnitude higher than for the cases considered in Section [3j here orbital decay occurs more than one 
hundred times faster during the high-eccentricity phase for a given value of Qb (note that for this example we have used a 
Q- value 10 times as large to clearly demonstrate Kozai cycles) . The effect of this strong tidal damping is to quickly reduce the 
mutual inclination to a value for which Kozai oscillations no longer occur (around 53°; see discussion below) and the apsidal 
lines become locked with 77 ~ 25-7T (panel (f)). The start of the slower phase corresponds to a value of planet b's semimajor 
axis of around 0.041 AU (the observed value is 0.0426 AU). Thereafter there is a slow decline in %, aj and ej,. 
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Figure 5. The Kozai regime. Here the initial mutual inclination and eccentricity of planet b are 60° and zero respectively, and Qi, = 100. 
The eccentricity of planet b is forced to high values initially (panel (b)) and there are corresponding variations in ij (panel (e)), however, 
this extreme behaviour comes at the price of severe tidal decay of the orbit of planet b (panel (c)) and the mutual inclination (panel (d)) 
until the system is no longer capable of driving Kozai oscillations. This occurs at around if, = 53° at which point the system becomes 



trapped on a limit cycle for which rj - 
summarized in Figure [6] Note that m c 



25tt 



and 



; m™ n /cos 60°. 



'/at 1 

mm 



0.02. This value of ij, = ip (see Table [3J is consistent with the analysis 



IN 

d 




Figure 6. (a): Maximum (Kozai) eccentricity when only quadrupole terms are included as a function of inclination and 7. The dashed 
black curve corresponds to 7 = 0, while the solid black curve, appropriate to the HAT-P-13 system with fe;, = 0.3, corresponds to systems 
with m c = m ™ n / 1 cos ij, | so that -f(m c = m™ m /\ cos = 7(m c = m™ In ) • | cos if, | (see equations l[2]l-(|4]l) decreases along both sides 
of the curve to zero at i;, = 90°. The maximum value of 7 is 4.2 at the end-points i;, = 54° and 126°. The circles correspond to the 
"true" maximum values of e;, for the fcj, = 0.3 case, calculated using the averaged code so that all tidal and octopole terms are included, 
(b): Orbit decay timescale as a function of if,. The solid curve and circles correspond to the solid curve and circles in panel (a) with 
Qi, = 10 5 , while the dashed and dotted curves corresponds to the same but with Qj, = 10 6 and Qi, = 10 7 respectively. We conclude from 
this that the HAT-P-13 system is unlikely to have a mutual inclination between 54° and 126° for reasonable values of Qi,. (c): The range 
of inclinations for which Kozai oscillations of ej, occur given a particular value of 7. 



Fabrvckv fc Tremainel (2007) give an expression for the maximum eccentricity, , attained during a Kozai cycle in the 



presence of other perturbing forces such as spin, tides and relativity. Their equation (34) (which holds in the case that the 
initial inner eccentricity is zero) can be rearranged so that 



ef = \/l-z 2 , (23) 
where x is the minimum positive root of the cubic 

x 3 + x 2 - (Ci + C 2 )x -Ci = 0, (24) 

with Ci = I cos 2 it, (0) and Ci — §7, i&(0) being the initial inclination of planet b (which equals the initial relative inclination 
since i c (0) = 0). Note that this expression does not include any contribution from octopole terms in the disturbing function, 
nor does it include terms of order e 2 . or higher in the tidal distortion. Figure |6ja) shows the dependence of e£ on it (0) (solid 
curve), taking into account the fact that the minimum mass of planet c is scaled by (coszt,) -1 so that the apsidal advance of 
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planet b is completely dominated by planet c when i c — 90° (see also Figure 3 of lFabrvckv fc Tremainel l |2007h ). The circles 



correspond to "true" maximum values of et for this case, calculated using the averaged code (with no dissipation) so that all 
tidal and octopole terms are included. While the maximum and minimum values of et do not vary from one Kozai cycle to the 
next when only point-mass quadrupole terms are included in the equations of motion, this is not true when extra accelerations 
are included, although the variation is generally small (less than 10%). We therefore integrated for several modulation cycles 
and recorded the maximum value of eb during that time. The results suggest that the terms not included in (|23|l make a 
significant difference for high inclinations. Also shown in panel (a) is when 7 = (dashed curve). Figure [6jh) shows the 
dependence of the orbit decay timescale on mutual inclination, for which an estimate is given by 

T a = (e I b i )~ 2 T circ . (25) 

That this gives a good estimate (rather than using, say, the average value of et during a Kozai cycle) is supported by the 
example in Figure[5]for which panel (c) provides an estimate over the first 10 4 yr of r a = 1.25 x 10 5 yr while (|25[) with e^ = 0.4 
(panel (b)) gives r a = 2.5 x 10 5 yr, that is, (I25[) represents an upper bound. The circles in panel (b) of Figure [6] correspond 
to those in panel (a) and show that (|23[) may be used in (|25[) to estimate r a in spite of the fact that not all relevant terms are 
included. We may therefore conclude from panel (b) that given an estimated age of around 5 Gyr, the mutual inclination of 
the HAT-P-13 system cannot be between 54 and 126 deg for values of Qb less than 10 6 . Note that if Qb is as high as 10 7 , the 
system will not yet have relaxed to the limit-cycle state. 

Finally, in order to be able to place bounds on the mutual inclination of any given observed system, it is important to 
know the range of inclinations for which Kozai oscillations of et occur given a particular value of 7. Figure Efc) shows that 
Kozai oscillations are completely suppressed for 7 >, 9 for any inclinations. This also has repercussions for the maximum 
possible stellar obliquity of a system, discussed in Section [6] 



3.3 Limit-cycle behaviour of retrograde systems 



While relaxed prograde systems are characterized by the libration of the angle r\ — tub — zo c , relaxed retrograde systems are 
characterized by libration of the angle ( — vot + vj c — 2Q.b- This can be understood as follows. 

Noting that r\ — 2ujb — (,, and that 77 ~ 2uib when £ librates, equations (f9|)- (|lip can be written 



e b = - [W T e b + Wo e c sin Q + A e b sin(2uj b t) + B r e c sin(2wbt), 
and 



with 



w: - w, 



°(t) 



cos C, 



+ A cos(2oj b t) + B r ( -^J cos(2oj b t), 



UJb — C — + Ob — —■UJc + ^b, 

where 



and 



7 / • \ ( mb\ a b -1 r /■ \ , a 
Vm c / V a c 



Wn 



W: = \(?£)e- 2 g2 {ib)Wn, 

with Wn defined in (I14p . Parameters which are zero for retrograde coplanar systems are 

5 ( a b " 



A 



while 



sin ib Wn and B 



cos ib + 




Wn- 



(26) 
(27) 
(28) 
(29) 
(30) 

(31) 
(32) 



Note that the rate of precession of the node of planet b about the invariable plane normal is negative for prograde orbits and 
positive for retrograde orbits (equation (|A15[0 . Given the symmetry properties of the various inclination functions (listed in the 
paragraph following (|A9|> ). we see that the relaxation behaviour of those retrograde systems for which £ librates is identical 
to that of prograde systems for which r\ librates (albeit with a slightly different limit cycle frequency). In particular, the 
theory developed in Sections 13. II and 13. 51 for prograde systems caries directly over to retrograde systems if one replaces f2(ib), 



fs(ib) and fi(ib) with g2(i b ), hz(ib) and hi(ib) respectively, as well as tub with —u)b- In particular, e b av \ib) = e b av> (ir — ib) and 



Aic(ib) = Aic(n— i b ), while (di b {ib)/dt) = 
because while fi(i b ) = /i(tt — ib), cos ib 



-{dib(n~i b ) jdt) (due to the factor sin(2i(,)). Note, however, that Cobiib) 7^ — Wi,(iT— ib) 
— cos(-7r — i b ). 
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Table 3. Data for cooling and Love number fitting laws. 
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Figure 7. Fitting functions for equations (133 I t and i l . i 4 1 ) for various values of the core mass of planet b. (a): Log-log relations between 
the tidal power required to maintain the equ ilibrium radius R^ eq ^ / Rj given the coolin g rate —Cb for such a structure, for core masses 
• 0, 40, 80 and 120 M®. Data is from Batygin, Bodcnhcimcr & Laughlin (2009) (circles) except for the endpoints corresponding 
l.OXRj which are chosen for best fit and so that the curves have the property that lim^-^Rj Cb = 0. (b): Exponential fitting 

->(<=?) 



(«!) 



to R 

laws for fc;, as a function of RY" 1 ' /Rj for the same selection of core masses as in (a). 



3.4 



The relationship between ei , the equilibrium radius of planet b, its Q-value and its core mass 



The above analysis implicitly assumes that the value taken for Rb is its equilibrium value. In fact, for a given set of orbital 
parameters and planetary core mass, the equilibrium value of Rb, R b eq \ depends on Qt,. This, in turn, determines ' for 
a given value of %. Moreover, the Love number depends on the core mass as well as the radius. We can quantify these 
relationships using the data provided in Table 1 of Batygin. B odcnhcim er fc Laug hlin ( 200 ^) a s follows 
For each value of the core mass of planet b, m CO re, 



(iZj**^ /Rj,Etide), where E tit i a is the tidal power required to maintain the radius at R K ^ ql given the cooling rate —Cb for such a 



Batygin. Bodenheimer fc Laughlin! (|2009h provide three sets of values of 



structure! The three values of fl- /*, correspond to the best-fit observed vame ±1., these being (1.20,1.28,1.36). Desirable 
properties of the relationship between Cb = Eud e and R b eq ^ /Rj are limn b ^Rj Cb = and d 2 (log 10 Cb)/d(log Rb) 2 < (see 
Figure 3 of iBodenheimer. Lin fc Mardlind (|2001h ). A function with these properties is 



j(e<?) 



\og 10 (Ri eq) /R.j] 



(33) 



with least-squares parameters b n and c n depending on m core , as well as the choice for Cb(R b eq ^ / Rj = 1.01)/Lq = A n (see 
Figure [71 a)). These are listed in Table [3] Also listed are least-squares exponential-fit parameters for the Love number for each 
core mass such that 



(34) 



Fitting laws (|33[) and (|34[) are plotted in Figure □ Following [Mardling fc Linl feooj ) Section 4, the rate of change of the radi 
of planet b is given by 



Rb 



— C b + Ef, 



-Cb/\E b \ +T- 1 

\Gm\l 'R b + a b m b Rlnl ~ (mb/m*)(a b /R b )' 



(35) 



4 Note that stellar insolation is also included in their calculations so our values for e! "' below are slightly overestimated 
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where at, is the moment of inertia coefficient of planet tjf] and Et = — \Gm r rribl a,b is the orbital binding energy, while to 
O(el), 0(f2e/«|) and 0(Slq/n%), the orbit-averaged expression for Eude is given bjjf] 



kb\ ( m*\ / Rb 
nib J \ ab 



9i\ 2 + (9l x 1 

rib I \rib 



+ (^] + 21e 2 



(36) 



Here f2 e = fib ■ ej, and fi 9 = fij, ■ qt, where Sib is the spin vector of planet b and q& = hf, x et, with et and hf, unit vectors 
in the direction of periastron and planet b's orbital angular momentum per unit mass respectively. Note that (|36|) assumes 
fib ■ hb — fib, that is, planet b is synchronously rotating with its orbit. The quantities fi e and fl q can remain non-zero for 
much longer than the tidal damping timescale if fl* x hj / and/or ht x h c / 0, where fl, is the star's spin vector and 
h c is the orbital angular momentum per unit mass of planet c. The terms involving fl e and fl q in (|36[) are associated with 



the so-called obliquity tide (jFabrvckv. Johnson fc Goo dman 2007); for the HAT-P-13 system they contribute less than 1% to 



{Eude} for all relative inclinations studied numerically, and will from hereon be ignored. The planet's radius will increase or 
decrease according to whether Eude is greater or less than Lb, at a rate enhanced by the factor (m,/mb)(Rb/a,b) which is 21.5 
for the HAT-P-13 system. Once the system achieves equilibrium (ie, once Eude ~ C-b), the radius of planet b will shrink on 
the timescale r , that is, the system will remain in the equilibrium state appropriate for the current value of ab- An example 
illustrating this behaviour is shown in panel (i) of Figure 1101 

For fixed Qb and m core , (|33[) and (|36[) may be equated to give a relationship between et, and i?^ 9 ' (using (|34[1 to express kb 



in terms of Rb)- Equation (|19[) provides a second relationship between these two variables for fixed ib, and combining the two 
gives one between R^^ and Qb- Solving this for R^^ for 10 2 ^ Qb ^ 10 6 and substituting the resulting values into II 191) allows 
us to plot e£ as a function of i?^ 9 ' and Qb for each value of m core /M® presented in Batygin. Bodenheimer fc Laughlin 



I 2009h . Such curves are plotted in panel (a) of Figure [8] for ib — and 30°, and panel (b) for ib = 40°, 45° and 50°. Also 



plotted (blue-dashed box) is the current la range of values of ' and R\^ q \ Since ib > ij ' for the curves in panel (b), 
the expression used for b's eccentricity is e^ a "' = ^(e™ n + e b rlax ), where e™ n and e b riax are taken from Table [2] Given 
T e tx T a oc e~ 2 (Section [2]) , we were able to confirm numerically that using e b av ^ rather than e™ a:c gives a good estimate for 
the decay timescale of ej,. This is in contrast to systems with high eccentricities (see Section \'i.2\ . 

Figure EJc) presents an alternative view for low inclinations, with each set of curves corresponding to a different core 
mass. Given the lower bound placed on Qb by r e (Section [2]), consist ent with t he fact that higher rather than lower values of 
Qb are suggested by the orbital parameters of short-period planets (Wu 20031 ) . Figure [8] suggests that the orbits of planets 



b and c are likely to be either near coplanar (mutually prograde or retrograde; see Section f3.3p . or have mutual inclinations 
between around 45° and 50° (or between 130° and 135°). However, a consistency argument due to Daniel Fabrycky (private 
communication) rules out the 45 — 50° and near retrograde coplaner cases. The steps in the argument are as follows, where 
we take the mutual inclination of orbits b and c, % c , to be 50° for definiteness: 

(i) LOb - oJe = Aoj = 4 ± 46 = -42 - 50°; 

(ii) The angle n = Ub — uj c + fib — fle librates around 180° for orbits with ib c = 50°, with a libration amplitude of around 
50°. Thus from (i), fl b - fie = Afl = 80 - 272° so that -1 sC cos Afl sC 0.17; 

(iii) Measuring inclinations with respect to the line of sight (as opposed to the invariable plane normal), ib = 83.4° ± 0.6; 

(iv) cos ibc = sin ib sini c cos(flb — fie) + cos ib cos i c ~ sin i c cos Afl; 

(v) Since ^ sini c ^ 1, cosii, c ^ cosAfi ^ 1 from (iv). But since cosibe = cos 50° = 0.64, this contradicts (ii). Thus a 
mutual inclination of 50° is inconsistent with observations. 

The argument against coplanar retrograde systems follows similarly, with cosi(, c ~ — 1 and the librating angle £ = u>b + u c — 
Afl = so that Afl = —48 - 44° with cos Afl > 0, making sin i c < 0. 

Figure [S] also suggests that lower rather than higher core masses are favoured. Note that for inclined systems for which 
ib < ij C1 ^, Te can be estimated by replacing Ao with A in JB), where A is given by (|12[) . Thus since A ^ Ao, the lower bound 
for Qb given by r e for coplanar systems represents a lower bound for higher inclinations in this range. 

More accurate measurements of et, and Rb will allow refinement of the statements above. Note that a small value of ib 
does not imply that a Rossiter-McLaughlin measurement of the stellar obliquity will result in a small value; its maximum 
depends not only on the mutual inclination of the two orbits, but also on the stellar obliquity relative to orbit c (see Sectionfi]). 



5 Formally ab is also a function of Rb, however, we take it to be constant at 0.26 in the numerical integrations presented in Section [4] 
corresponding to the moment of inertia coefficient of a polytrope of index 1. 

6 Note that equation (73) in lMardling fc Lin (2002) should read {Etide) n °t {Etot}- Note also that their ki refers to the apsidal motion 
constant which is half the Love number. 
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Figure 8. The relationship between e b av \ the equilibrium radius of planet b, its Q- value and its core mass for the HAT-P-13 system 
(note different scales on panel (c)). (a): Sets of curves for if, = (black) and if, = 30° (red), with the core mass increasing from bottom 
(nicore = 0) to top (m core = 120M®), and 7713 = m™ %n / cos if, . The open circles correspond from right to left to Qb = 10 2 , 10 3 , 10 4 , 
10 5 and 10 6 ; notice how they tend to bunch up as Qb increases and Ry 9 ' — > 1. The blue dashed box indicates the current la range 
of values of e[°"' and Hj . (b): Sets of curves for if, = 40 (black), if, = 45° (red) and if, = 50° (green dot-dashed), with core masses 

and Q- values the same as in panel (a). Since if, > the expression used for the eccentricity of planet b is = \{&™ m + £™ ax ), 

where e™ m and e b nax are taken from Tabled Given the lower bound placed on Qj, by r e (Section |2J, panels (a) and (b) suggest that 
for mutually prograde systems, the orbits of planets b and c are likely to be near coplanar, or have mutual inclinations between around 
45° and 50°. Lower rather than higher core masses are also favoured, (c): A alternative view. Sets of curves for m core = OMj (black), 
4OA/0 (red dashed) and 80 Mq (blue) for (from top to bottom of each set) 0°, 10° and 15°, and (circles from right to left along each 
curve) Q b = 10 4 , 10 5 , 10 6 and 10 7 . The green dashed box indicates current la range of values of e^""' and R b eq \ 



3.5 Timescale for decay of the mutual inclination in a relaxed system 

Taking if, as a proxy for the mutual inclination (since we are assuming most of the angular momentum of the system resides 
in the outer orbit), we can use (|A13f) to estimate a timescale for the decay of the mutual inclination of a relaxed system. 
Figure |9j a) shows 10 years of evolution of if, for the systems shown in Figure [2] for which Qb = 10 and kb = 0.3. In order to 
clearly demonstrate long-term trends, we have plotted the quantity if, — it, c (0), where if, is measured relative to the invariable 
plane and ib c (0) is the initial mutual inclination. In particular, note that the general trend is positive for if, = 10° and 20°, 
while for if, = 30°, if, = 40° and 50° it is negative. For systems for which 77 librates (here, if, = 10°, 20° and 30°), this can be 
understood as follows. 

Consider equation (|A13[) for the rate of change of if,. The current slope of the trend may be determined by taking a time 
average of dib/dt holding e c , if, and ab/a c constant. Replacing ej and sin n by their limit-cycle counterparts e£ and r/^' 
respectively (equation (|17[) ). and recalling that Wb + ro c — 2fif, = 2u>b — Tj, we obtain 

/ di b \ _ 1 f T di b 15 / m c \ / a b \ 4 _ 5 r av ) r av) . r. ,. , l**i.f\\ 

\-dtl = t^LtJ -tt=~32 nb {^n-){-J £c V ^^[AW-^hiW] 
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Figure 9. (a) Decay (or increase) of ib for Q = 10 and fcj, = 0.3, where ij, is measured relative to the invariable plane and i(, c (0) is 
the initial mutual inclination. Solid black curves represent numerical integrations, red dashed curves give the trend according to l |37|l 
and green dotted curves are horizontal reference curves to guide the eye. Trend curves are not given for i(, c (0) = 40° and 50°, systems 
for which the theory is no longer valid. Note that for the cases if, c (0) = 40° and 50° the oscillatory behaviour is dominated by the 
quadrupolc contribution to di^/dt (which is proportional to sini c and has frequency — S7 C ) while for it{0) = 10°, 20° and 30° the 
octopole contribution is evident (see equation l)A13|l 'l. (b) The function *£(ib', A, 75^) defined in equation l |38|l . Note that ^ > for % < ijj, 
consistent with numerical solutions. 



25 ( a,b\ 2 2 -4, Tf /. a 3\ -1 
= 32 \~) 6c£c *W A >T e J ■'Tare 

= (ib-iDr' 1 , (37) 

where 

A, 7 Sc) = /a(i»)an(2i») [/ 4 (*0 - ^At c (i b )h 4 (ic)] A" 2 (38) 

with A defined in (I12p . i£ is the first root of and Tj is time it takes for the system to relax to % = i£. For ij, = 30°, 
Ti — 80 Gyr. The function ^(if,; A, 7ejj) is plotted in Figure [9jb) as a function of it, for HAT-P-13 system parameters with 
kb = 0.3. We see that $ > for 4 < i% ~ 24°, that is, if, actually increases for this range of inclinations, while for greater 
inclinations ib decreases. The red dashed lines in Figure |9ja) have slopes equal to (dib/dt), and are in good agreement with 
the general trend of the numerical solutions. Also shown are numerical solutions for ib = 40° and 50°. The slope of the trend 
for i b = 40° is 0.4(10 5 /Qi>)° Gyr -1 , while that for i b = 50° is a factor of three higher. Thus the timescale for the decay of 
the mutual inclination of HAT-P-13-like systems for which if, < 50° is considerably longer than the age of the system for 
reasonable values of Q b . 

A similar analysis can be done for retrograde systems, and we can conclude generally that inclined systems cannot relax 
to the coplanar prograde or retrograde state as long as e b av ^ > 0, and that they relax to a mutual inclination given by one 
of the roots of \&(ib) = (or 9(tt — if,) = for retrograde) as long as r, < r c and r, < r a . Note that the mutual inclination 
changes by the same amount as %, and that (|37[) does not apply to systems with mutual inclinations greater than i^ C1 ' for 
which it represents a lower bound. 



4 A SCENARIO FOR THE ORIGIN OF HAT-P-13-LIKE SYSTEMS 

The high eccentricity of HAT-P-13c suggests either a violent scattering history, or Kozai-type interactions with another planet 
or star, or that it was formed through gravitational c ollapse; it seems unl ikely that such a high eccentricity could result from 



single planet-disk interactions alone (see, for example, Artvmowicz ( 19931 )'). Kozai forcing would also affect planet b, and it is 



unclear without a detailed study whether or not a suitable configuration exists which would not cause the rapid decay of its 
orbit. 

Here we focus on the scattering scenario, introducing a third planet ("planet d") which is ultimately ejected from the 
system. The mass ratio of planet b to planet c is too low for c to have attained its high eccentricity during a scattering event 
with b. Three models are presented, all of which have identical initial conditions except for the eccentricity of planet d, ea, 
which for model 2 differs from 1 by two parts in 10 J while for model 3 it differs by 0.05. The outcomes are significantly 
different, with model 1 producing a value for e c almost identical to that for HAT-P-13c, while models 2 and 3 produce lower 
values at 0.46 and 0.36. Of particular interest is the relative inclination of the orbits of b and c following the escape of d, and 
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Table 4. Data for models 1, 2 and 3 at t = and around time planet b first achieves equilibrium 



Qb = 40 



model 1 



model 2 



model 3 



* (yr) 





1.5 Myr 





10 Myr 





35 Myr^ 


m b (Mj) 
m c (Mj) 
m d (Mj) 


0.85 
15.2 
12 




0.85 
15.2 
12 




0.85 
15.2 
12 




a b (AU) 
a c (AU) 
a d (AU) 
Rb (Rj) 


0.043 
1.7 
3.097 
1.5 


0.0406 
1.17 

oo 

2.34 


0.043 
1.7 
3.097 
1.5 


0.0422 
1.16 

oo 

1.44 


0.043 
1.7 
3.097 
1.5 


0.04255 
1.17 

oo 

1.14 


eb 
e c 


0.005 
0.1 


0.007 
0.68 

1 .U-LO 


0.005 
0.1 

n 0^009 

U.UJUUZ 


0.004 

0. 45 

1. UJ 


0.005 
0.1 
n i 

U. 1 


0.001 
0.36 

1 m 


tp*b 

'4>,c ^ 8, 


0° 
15° 


3.5 - 18° 
7-8° 


0° 
15° 


3 - 40° 
18 - 21° 


0° 
15° 


23 - 53° 
12 - 19° 


h 


10° 


6 - 14.5° 


10° 


8 - 32° 


10° 


30 - 38° 


i« 


5° 


4° 


5° 


11.3° 


5° 


4° 


*f 


5° 


2° 


5° 


6.2° 


5° 


2.3° 


"fa DC 


15° 


10.2 - 10.6° 


15° 


19.5 - 21° 


15° 


34 - 35.5° 


dip 




0.2° 

4° 




0.2° 

11.2° 




0.35° 

3.8° 



a Equilibrium has not yet been established. 



the accompanying stellar obliquity relative to planet b, both of which differ significantly from model to model, as does the 
equilibrium value of Rb- 

The initial configuration data are listed in Table 3] together with those at the time planet b first achieves an equilibrium 
radius, that is an equilibrium between the rate at which tidal energy is injected and the rate at which the planet can cool. 
Data of particular significance are highlighted in bold, fnclinations with superscripts (i) and (/) are measured with respect to 
the original and final invariable planes respectively, the former including planet d and the latter not. The qu antity 8ip is the 



angle between the invariable planes before and after the escape of planet d. Again we use the averaged code of lMardling fc Lin 
(2002), ;his time without suppressing the evolution of the radius of planet b which is evolved according to the scheme discussed 
here in Section [3.41 with m CO re = 80M®. Note that this code involves averaging over the innermost orbit only; the two outer 
orbits are integr ated directly and are therefore susceptible to instability through the overlap of mean-motion resonances 
I Mardlinj l2009al ) . While the mass of planet b is taken as that of HAT-P-13b, the mass of planet c is taken as the minimum 
mass of HAT-P-13c, that is, it is not scaled by cosij,. The mass of planet d is taken to be similar but less than that of planet c, 
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ensuring that the probability that planet c is left with a high eccentricity is significant. Initial values of the eccentricities and 
inclinations are chosen to be consistent with formation and subsequent migration in a protoplanetary disk, and the period 
ratio (which is similar to that for Jupiter and Saturn) puts th e system just ins ide the stability boundary of a system with 



these masses and eccentricities (see, for example, Figure 9(a) of lMardlinel |2009ah ). consistent with having just emerged from 



the protection of the diskQ Moreover, the initial values of a c and the semimajor axis of planet d, ad, are such that the energy 
needed for escape of planet d and provided by the orbit of planet c reduces a c to a value similar to that of HAT-P-13c. The 
initial radiu s of planet b is 1.5 R j, consistent wit h the radius of a young planet recently arrived at its present location (see 



Figure 1 of Bodenheimer. Lin fc Mardlin j ( 2001 ); disk lifetimes suggest that such a planet would probably have a higher 
radius). The initial mean longitudes and orientation angles of planets b and c are specified relative to the outer orbit; these 
are Xbd ~ 0°, ibd = 5°, uibd = 0°, Qbd — 0° and A c< j = 180°, i c d ~ 10°, oj c d = 0°, Q c d = 180° respectively and are such that 
longitudes are measured with respect to the periastron direction of orbit d. 

The time evolution of various quantities is plotted in Figures 1101 1111 and 1121 with detailed descriptions and discussion 
provided in the captions. The choice of an initial period ratio for orbits c and d of around 5:2 affects possible outcomes in the 
following ways (given that we wish to approximately reproduce the HAT-P-13 system). In view of the fact that we wish to 
place the system near the edge of the stability boundary, the choice of e<j(0) is restricted by the fact that a value significantly 
higher would render the system far from the boundary and hence violently unstable. Our choice of mj affects the choice of 
dc(O); a smaller value of rrid requires a smaller value of a c (0), however, it tends to produce smaller values of e c following 
escape of d and is less likely to produce high stellar obliquities (see model 3). On the other hand, a heavier companion is more 
likely to eject planet c (an outcome also possible for mj = \1Mj\ here we restrict ourselves to systems in which the outer 
planet is ejected). A general detailed study is required to qua ntify possible outcomes of such a scatte ring scenario, especially 
in the light of the recent discovery of two retrograde systems ( Winn et al. 20091 ; Anderson et al. 20101 ). For now our aim is to 



demonstrate its potential to produce a range of eccentricities, mutual inclinations, stellar obliquities and planetary radii. 

Of particular interest is the stell ar obliquity relative to th e orbit of planet b, ?/>,(,, a quantity which can be measured 
directly (at least in sky projection; see Fabrvckv fc Winn ( 20091 ) for a discussion of the statistical properties of this quantity). 



In the following section we outline the mechanics of stellar obliquity in a two-planet system. 
4.1 Stellar obliquity in a two-planet system 

The stellar obliquity relative to each planetary orbit in a system reflects conditions at the time of its formation. In the 
scattering models presented in this section we assume zero stellar obliquity relative to planet b initially, while the stellar 
obliquity relative to planet c is 15°. Table [4] lists the ranges for xj)*b following the escape of planet d, demonstrating that 
significantly different outcomes are possible from models with very similar initial conditions. The variation in tp^j, depends 
on two quantities: the angle between the stellar spin axis and the normal to the invariable plane, 9,, and the angle between 
planet b's orbit normal and the normal to the invariable plane, %. When i/>„f, 7^ 0, the variable torque on the spin bulge of the 
star from planet b results in nutation of the star's spin axis, that is, a variation of 9*. This can be quite significant; in model 3 
it is 7° compared with 1° for model 1 (see Table [4])[f] Its average value, however, depends on the history of the system; if no 
planets have been ejected since the system's formation, 9, is likely to be modest, while the escape of one or more planets can, 
depending on how much angular momentum the planets carry away, result in an i nvariable plane normal which points in a 
significantly different direction to the original (Oip in Table [4|, thereby affecting Barker fc Ogilvie ( 20091 ) have shown that 



the decay timescale for ip*b is around T a , the timescale for the decay of the orbit, so very little reduction in the average value 
of 9, is expected over the lifetime of the orbit of HAT-P-13b. 

Regarding the relevant angles as spherical polar angles, ip^i is given in terms of 9*, the spin axis node angle, as well 
as ib and Qb, by 

costp*b = ^* ■ hj, = sin#„ sin % cos((p* — Qb) + cos#* cos ib- (39) 

Since the variations in 9* and % are small (but see next section), tji^j, cycles approximately between \ib + 9,\ and \ib — 9,\ 
as </3* — Qb cycles between and 2-7T over b's precession cycled or, since the orbit of planet c contains most of the angular 
momentum of the system so that 9, ~ ^>* c , the variation is approximately \ib ± ip* c \- In fact, since the nutation period of the 
star is equal to the precession period of planet b, extrema of ip*b coincide with extrema of if} tc (see panel (e) of Figures [lOl to 
E]). 

A significant difference between the first two models and model 3 is that before the escape of planet d, the latter suffers 
a significant transfer of angular momentum between orbits c and b during a period of high eccentricity of planet c, resulting 

7 The disk will still be present at this stage, but its density will not be high enough to suppress eccentricity growth and prevent the 
overlap of neighbouring mean-motion resonances. 

8 Note that the moment of inertia of the star is about three times that of b's orbit. 

9 Note that the precession rate of the spin axis node is more than 15 times slower than that of b's orbit for the examples in this section. 
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Figure 10. A chaotic origin for the HAT-P-13 system (note different timescales for each panel). At time t = three planets are present, 
with migration leaving planet b 0.043 AU from the star with a radius of 1.5-Rj and a core of mass 80Afgj, and planets c and d with period 
ratio 2.46. While the presence of an outer disk has previously limited the variation of the eccentricities thereby protecting the system 
against instability, by t = the disk surface density has reduced sufficiently to allow the system to become unstable. Panel (a) shows all 
three eccentricities, each remaining moderate for the first 6000 years. After 16,000 years planet d escapes the system. Following escape, 
the orbit of planet b precesses around the new invariable plane whose normal is approximately parallel to that of planet c (since with 
d removed from the system it now contains 99% of the total angular momentum; compare panels (b) and (c)). The inclinations of the 
three orbits to the original invariable plane are modest initially, and they remain so during and after the scattering process. Following 
the escape of planet d, the stellar obliquity relative to b's orbit (panel (e), black solid curve) oscillates about a mean equal to ib (blue 
dashed curve) and with an amplitude equal to ip* c , the stellar obliquity relative to c's orbit (red dashed curves; see also panel (f) and 
Figure fT3)l . Panel (g) shows the early behaviour of ej, ; before the scatter planet b is effectively decoupled from the rest of the system and its 
eccentricity monotonically decreases on the tidal circularization timescale, while after the scatter it is governed by planet c. An artificially 
low Q-value of 40 is used for planet b, and its radius and Love number are evolved according to (1331) and I I 341) respectively (fc(,(0) = 0.13 
for -Rb(O) = 1.5-Rj and m CO re = 80Mq). Panels (h) and (i) show the long-term evolution of e;, and Rb respectively. The system evolves 
to a limit cycle after about 0.12 Myr, with e^ a "' decreasing as Rb increases (given the dependence of e^ v ^ on Rb via jt' de oc (Rb/at,) 5 ) 
until an equilibrium value of Rj>, R^^ = 2.33/2 /, is reached after 1.6 Myr. This value of R^^ corresponds to = 0.007 and agrees 
favorably with the theoretical values R^ q ^ = 2.12Rj and ej, = 0.008 obtained using the procedure described in Section 13.41 Note that 

more realistic (ie, higher) values for Qb would result in higher equilibrium values of e^ and lower values of R b eq ^ (see Figure [8}. During 
the 3.5 Myr integration, the semimajor axis decreases from 0.04300 to 0.04055 AU, giving an orbital decay timescale r a = a(,/d(, of 
158(Q(,/10 J ) Gyr. Also shown in panel (i) are curves Rb(t) = i?i,(0)exp(f/r-|_) (red dashed curve) and Rb(t) = -R(,(t* )exp[— (f — t*)/T a ] 
(blue dotted curve), where t+ = (mj, /m» ) (a;, / Rb)r a and t* = 3.5 Myr. These curves demonstrate that initially the rate of change of the 
radius is dominated by tidal heating, that is, the cooling term contributes very little (see equation 1 1351) ). and once equilibrium is reached 
it is continually reestablished as ab decreases on the timescale T a . 



in a particularly high maximum value of if},b- The importance of this difference can be understood as follows. In each model, 
the (average) value of ip* c depends on its value before the scatter, as well as the change in the inclination of c's orbit relative 
to the original invariable plane. Since the orbit of c effectively coincides with the invariable plane once d has left the system, 
ib — ip* c for a system in which there is no close encounter between b and c (at least for the cases considered here for which 
ip*b(0) = 0), and the variation in ip*b is approximately equal to 2tp* c - On the other hand, ib is significantly greater than ip* c 
in model 3 because of the close encounter of b and c, and as such the maximum value of ij)*b is high. 

Figure [13] illustrates the mechanics of stellar obliquity in a two-planet system for the cases (a) if>* c = 20° and ib = 70° (a 
configuration ruled out for the HAT-P-13 system because 54 < ib < 126°; see Section [3.2[) and (b) i[>* c = 70° and ib = 20°, 



Tidal evolution of inclined systems 19 



eccentricities 



inclinations re I to invariable plane 



inclinations rel to NEW invariable plane 



1 1 1 1 1 1 1 1 
- (□) 







y c 




VW b 



5000 1(T 1.5x1CT 2x10^ 
years 

semlmajor axis of planet c 




(c) 



stellar obliquity rel to orbit b 



1CT 2x10^ 3x10^ 
years 

stellar obliquity rel to orbit c 



III 


o 

cn 

CD 




1 1 1 1 1 








- il 








_D O 






1. I 1 .1. I 1 I 1 1 I 1 I 1 1 I 1 1 1 1 




m 


, , > , ,-- 





(f) 



5000 10' 1.5x10' 2x10' 

years 



10' 2x10' 3x10' 4x10' 

years 

detail of e b (long term) 




10' 2x10' 3x10' 4x10' 
years 

radius of planet b 






5000 10' 1.5x10' 2x10' 

years 

Figure 11. A different outcome for the HAT-P-13 system (note different timescales for each panel). Initial conditions are the same as 
for model 1 but with a difference in e j of 2 X 10 — 5 . Following the scatter of planet d from the system, e c = 0.45 (panel (a)), resulting in 
a smaller value of e b av ^ (panel (g)). This time the mutual inclination reaches 20°, as does the stellar obliquity relative to c's orbit (panels 
(c) and (f)). This results in a maximum stellar obliquity relative to b's orbit of i b + tptc = 40° and a minimum of 3° > i b — V'*c (panel 
(e); the variation of ip* b is not a simple sinusoid when its mimimum is near zero, similar to the behaviour of eccentricity in the same 
circumstances). Note that the precession timescale is longer than that of model 1 by a factor (e| /£c ) 3 , where the superscripts refers to 
the model numbers (see equation I I14II ). The long-term evolution of e b (panel (g)) is a limit cycle for which e; a "' increases slightly for the 
first 10 Myr or so, corresponding to the relatively rapid (but still very slow) decrease of R, b (panel (h)), then decreases once equilibrium 
is established. The theoretical prediction according to Section 13.41 gives R^ eq) = 1A2Rj and e[ av) = 0.0037 when a b = 0.0426. During 
the 20 Myr integration, the semimajor axis decreases from 0.0430 to 0.0422 AU, giving an orbital decay timescale r a = 2590(Qf,/10 ) 
Gyr. Also shown in panel (i) are curves R b (t) = -Rf,(0)exp(— t/r—) (red dashed curve) and R b (t) = i?(,(t*)exp[— (t — t t )/r a ] (blue dotted 
curve), where t_ = (m b /m t )(a b / R b )\E b \/ C b (R b (0)) and t, = 20 Myr. These curves demonstrate that the rate of change of the radius 
is never dominated either by heating or cooling (see equation 1351) . being close to equilibrium initially. After around 10 Myr equilibrium 
is established and R b decreases on the timescale r a . 



so that 50° ^ ip*b ^ 90° f° r both. Also indicated is the original invariable plane in the case that the system contained a 
third planet which has since escaped, and whose angular momentum was such that the stellar spin axis was parallel to the 
invariable plane normal. The ramifications for the origin of retrograde systems are clear (although configuration (b) would 
probably require a close encounter with a passing star rather than the escape of a companion planet to cause such a dramatic 
change in c's orbit). Given the amount of angular momentum transferred from c to b depends on all the system parameters, 
it is easily conceivable that a scattering scenario similar to the one described here (either bound or flyby, perhaps involving 
exchange of c and d) could produce retrograde systems such as HAT-P-7 and WASP-17. It is interesting to note that it is 
possible for the relative inclination of a system to pass through 90° without destroying planet b through Kozai oscillations as 
long as 7 is large enough (see Section f3. 2[) . 

We end by considering the effect on the mutual inclination of non-zero stellar obliquity with respect to the invariable 
plane, which can be significant if 9, is significantly non-zero. 
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Figure 12. High stellar obliquity due to a close encounter between b and c (note different timescales for each panel). Initial conditions 
are the same as for model 1 but with = 0.1. This time the mutual inclination reaches 35°, while the maximum stellar obliquity relative 
to c's orbit after the scatter is 19° (panels (c) and (f)). This results in a maximum stellar obliquity relative to b's orbit of if, + ip* c = 53° 
and a minimum of if, — ip* c = 23° (panel (e)). Unlike models 1 and 2 where the stellar obliquity can be attributed almost entirely to 
angular momentum transfer between c and d (and the initial value of if)* c ), if, and V*c are significantly different, a result of a close 
encounter between planets b and c around 17,000 yr when e c reaches 0.87 (see panel (g) in which the periastron separation, p c , reaches 
a minimum of 4a;,). During this high-eccentricity phase, the torque from planet c tilts the orbit of planet b through 30°, and while some 
of the associated angular momentum is subsequently returned to c's orbit, the mutual inclination remains high following the escape of 
planet d. Since the value of e c is relatively low at 0.36 after the escape (panel (a)), e£ is a mere 0.001 (panel (h)), with a similarly 
small amplitude in spite of the relatively high mutual inclination. As a consequence, tidal heating is weak and the predicted equilibrium 
radius of planet b is only 1.14/?/ (recall Qf, = 40). Panel (i) shows the evolution of Rf, together with curves Rb{t) = i?.f,(0)exp(— t/r_) 
(red dashed curve) and Rb(t) = /?{,(£* )exp[— (t — t*)/T a ] (blue dotted curve), where this time t, = 35 Myr. These curves demonstrate 
that initially the rate of change of the radius is dominated by cooling (see equation 1351) . and does not appear to reach equilibrium during 
the time shown. 



4-. 1.1 Effect on the mutual inclination of non-zero stellar obliquity with respect to the invariable plane 



The numerical results presented in Section 3 assume zero stellar obliquity relative to the invariable plane, that is, 9, = 0. 
Consequences of this are that the relaxed state illustrated in Figure[T](red curves) and Figure[2ja) exhibits constant amplitude 
variations, and that the inclination if, is effectively constant. The latter is clearly not the case in the three models discussed 
above (panel (c) of Figures [TOl 1111 and I12|) . Figure [I4f a) shows how the amplitudes of variation of and if, are modulated 
when 9* is non-zero. Here 9* — 50°, if, = 30° and Qt = 10, and the initial values of ef, and n are ej""' and zero respectively. 
The limit-cycle amplitude varies noticeably and the modulation period is slightly longer, with the behaviour persisting during 
a 10 6 year integration. This can be understood as follows. The torque on the orbit of planet b due to the star's spin oblateness 
is given by equation (48) in lMardling fc Lin (2002), and this contributes to the rate of change of it according t o equation (29) 



of the same paper. Expressing the stellar spin vector J2« and the basis vectors ej, qt and hf, referred to in Mardling fc Lin 
(2002 m terms of the invariable plane reference basis via the star's node angle (p», and the Euler angles uib, ^f, and if,, 



the contribution to dib/dt from the star's spin oblateness becomes 



See also Section 13.41 here. 
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Figure 13. Illustration of high stellar obliquity in cases where (a) planet c is scattered through a small angle during its interaction with 
a third planet, and planet b is scattered through a large angle during its interaction with planet c, and (b) planet c is scattered through 
a large angle during its interaction with a third planet (or passing star), and planet b is also scattered through a large angle during its 
interaction with planet c, but in such a way that it ends up with a small value of if,. Here we are assuming that the stellar spin direction 
is parallel to the invariable plane normal before interaction with the fourth body, and that the system was originally coplanar. While 
both configurations have the same range of values of ip*b = |*bi^*c|i configuration (a) is ruled out for HAT-P-13 because 54 < if, < 126° 
(see Section |3.2|I . while configuration (b) cannot be ruled out (although it would probably require a close encounter with a passing star 
rather than the escape of a companion planet to cause such a dramatic change in c's orbit). 



LO 




2x10 4 4x10 4 ™ 2x10 4 4x10 4 

time (yr) time (yr) 



Figure 14. The effect of non-zero stellar obliquity relative to the invariable plane on the relaxed state for (a) ef, and (b) if,. The 
limit-cycle amplitude (ie, the amplitude of variation of ef,) varies regularly on a timescale of twice the limit-cycle period, the latter being 
27r/2dj(,, and its modulation period is slightly longer. The amplitude of variation of if, is significantly enhanced by an amount proportional 
to sin(20*). For comparison, the variation of if, for the case 0*(O) = is shown (red curve in panel (b)). The small variation is due to 
nutation of the stellar spin axis through about 2.5°. 



— = Vbk* ( — -) ( — -] [cos #« cos it cos(f2{, — ip*) + ~ sin 8, sin if, sin(2fi[, — 2yj*)l . (40) 

at * 2 V a,b J \ rib ) J 

Taking (9*, if, and tp* to be approximately constant over an orbit precession cycle (the precession rate of the stellar spin axis 
is 40 times slower than the orbit precession rate for this example), and putting Qf, = Clf,t + fif,(0) where Clb is given by (|A15|I 
and is also approximately constant, the variation in it, due to the star's spin oblateness, Aib, is obtained by integrating (|40|l 
over an orbit precession cycle to give 

Aifc = 7 r n 4 S in(2V>*c) (41) 

where 7» pln is the contribution to 7 from the star's spin quadrupole moment and is given in Table Q] for HAT-P-13 for 
m c = m™ n . Putting yj tc — 50° then gives Alt, = 1.1°, in good agreement with Figure I13f b). Such a variation in turn 
produces a variation in e ™ ax < mm G f t ne order seen in panel (a) of the same figure. 



5 SUMMARY 

The ideas and results presented in this paper can be summarized as follows: 



22 Rosemary A. Mardling 



1. Generalizing the results of iMardlinj 1 2007 ) to non-coplanar systems for which most of the angular momentum of the 



system resides in the outer orbit, we find that under the action of tidal dissipation in planet b, the system evolves to a limit 
cycle rather than a fixed point in ef, — r\ space, with the average value of ej, e^ v \ decreasing and the limit cycle amplitude 
increasing with increasing mutual inclination, and limi 6 ->o e£ = e b eq \ For the HAT-P-13 system, limit cycle behaviour 
occurs for if, < 33° (libration of r\ around 2nn) and 46° < ib 54° (libration of r\ around (2n + No limit cycle exists for 
33° < ib ^ 46° (rj circulates). 

2. For systems with 54° < if, < 90° (90° < ib < 126°), Kozai oscillations coupled with tidal dissipation in planet b act to 
reduce (increase) the mutual inclination until ib < 54° (ib > 126°) on a timescale much less than the age of the system. We 
conclude that the HAT-P-13 system cannot have a mutual inclination between 54° and 126° for Qt £ 10 6 . 

3. For retrograde systems, the limit-cycle behaviour is the mirror image of that for prograde systems, apart from a slightly 
different limit-cycle frequency (compare (|32[) with (|16[0 . 



4. The analysis and conclusions of iBatygin. Bodenheimer fc Laughlin (2009) for the HAT-P-13 system are valid as long as 



the mutual inclination of planets b and c is less than around 10° (which may well be the case; see next point). For higher 
prograde values of if,, a measurement of ef, does not unambiguously determine kb, although one can make arguments about 
the likelihood of a system being near the top or bottom of the modulation cycle of ef,. Similar statements hold for mutually 
retrograde systems. 

5. We have derived a relationship between the average eccentricity, e b \ the equilibrium radius of planet b, R b eq \ its Q-value 
and its core mass (Figure [8]), and conclude that for Q- values greater than the lower bound imposed by the timescale for decay 
of e b av \ the orbits of planets b and c are likely to be either near prograde coplanar, or have mutual inclinations between 
around 130° and 135°. Lower rather than higher core masses are favoured. More accurate measurements of et and Rb will 
allow refinement of these statements. 

6. Inclined systems cannot relax to the coplanar prograde or retrograde state as long as e£ > 0, and instead relax to a 
mutual inclination given by one of the roots of ^(ib) = or ^(tt — if,) = (see equation (|38|) V This will occur as long as 
Ti < t c and Tj < r a , both true for the HAT-P-13 system, however, in this case, Ti is much greater than the age of the system. 

7. A viable formation scenario for the HAT-P-13 system is that it originally contained a third planet which was scattered out 
of the system when the protoplanetary disk density dropped below the critical level for stability. Such a scenario is capable of 
producing a high eccentricity for planet c as long as the mass of planet c is sufficiently high, and may also produce significant 
mutual inclination, stellar obliquity and inflated planetary radii. 

8. A Rossiter-McLaughlin measurement of the sky-projected stellar obliquity relative to planet b, together with a measurement 
of the mutual inclination, will allow us to constrain the stellar obliquity relative to planet c and hence obtain knowledge about 
the formation history. 

In conclusion, it is likely that many more HAT-P-13-like systems will be discovered in the future including systems in 
which the outer body is a binary star companion. Such systems will contribute a wealth of information not only about the 
internal structure of the short-period planet, but also about the formation history of the system. 

Note added in proof: 

Since this paper was submitted more refined data for the HAT-P-13 system have become available ( Winn et al. 2010T ). in 



particular, a new estimate for the eccentricity of planet b of 0.01421q'qo44- Moreover, there is now evidence for a distant third 
body in the system, as well as a Rossiter-McLaughlin measurement of the sky-projected stellar obliquity, the latter strongly 
suggesting that the stellar spin and the orbit normal of planet b are aligned. The conclusions drawn here remain valid, in 
particular those regarding coplanarity or otherwise when the refined value of co b — u) c is used. 
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APPENDIX A: SECULAR EQUATIONS FOR A NEWTONIAN POINT-MASS SYSTEM 

In Mardlin g (2009b), the secular equations governing the evolution of the orbital elements of the inner and outer binaries of 
a hierarchical triple are given to arbitrary order. These are derived using a spherical harmonic expansion of the disturbing 
function 1Z expressed in Jacobi coordinates, the latter having the dimensions of energy and denned to be such that the total 
energy is 

1 Gra.mi, lG(m,+m b )m c 

hj = It. (Al) 

Z ab 1 a c 

Noting from the numerical solution for the HAT-P-13 system presented in Section [3] that e b , ab/a c and max(sini c ) are all of 
order 0.01, and taking the invariable plane to be the reference plane, the orbit-averaged disturbing function to order sufficient 
in et, a b /a c and sini c to produce the dominant terms of each of the rates of change of the elements is 

■R. = n q +TZ , (A2) 
where the quadrupole and octopole contributions are given by 

TZ q = \nbdlnl (J^J (^j £c 3 {(1 + |et)/i(«b)/i(«c) + x e & sin 2 4 cos^o;;,) 

+ | cos(Oi, — fi c ) sin 2i b sin 2i c + | cos(2f2(, — 2f2 c ) sin 2 it sin 2 i c } + 0(x 5 ) (A3) 

and 

TZ = -i§/i b a 2 n! (^-^j (~~) e b e c e~ 5 {/ 2 (ii,) cos(tu b - zu c ) + Q2{ib) cos (zu b + zu c - 20b)} + 0(x 7 ), (A4) 



24 Rosemary A. Mardling 



respectively, where x represents any one ol et,, a b /a c or sini c (for example, x 2 may represent e b sini c , a nd cos i c = 1 + Q( x 2 )). 
Note that here we have assumed that mj/m, <C 1 and m c /m* <C 1 (in contrast, the equations in Mardlinel 1 2009b ) are 
completely general). The inclination functions fi(ib), fziib) and g 2 (h) are listed as follows together with fz(ib), fi(ib) and 
/i4(ifj) which appear in the evolution equations below, and hs(ib) which appears in the theory of relaxed retrograde orbits 

(A5) 
(A6) 
(A7) 
(A8) 
(A9) 
(A10) 
(All) 

The inclination functions (|A5[l - (|A9[) have the symmetry properties /i(it,) = /i (7T — it), fz(ib) = gzin — ib), fz{ib) = fe(7T — ib), 
fi(ib) = hiiji — ib), /n(0) = 1, 32(71") = h-j,(Ti) = hiiji) — 1 and 32(0) = /2(vr) = 0, all of which are relevant for the comparison 
of prograde and retrograde systems (see Section f3 . 3 [I . 

Given our aim of understanding the long-term behaviour of HAT-P-13-like systems and having demonstrated empirically 
that they maintain small values of sini c and evolve towards small values of e t , on a timescale equal to three times the tidal 
circularization timescale for any initial relative inclination, we follow iMardlina 1120071 ) and write down the equations governing 
the secular evolution of t he orbital elements, retaining only leading order terms in e&, ab/a c and sini c . Using (1A2[) . Lagrange's 
planetary equations are (jMurrav fc Dermottllioool l 



(Section EO)): 






fifa) = 


|(3 cos 2 ib — 1), 






.Hit) = 


g(l + cos ib) (15 cos 2 ib — 


10 cos if, 


-1) 


h{ib) = 


| (5 cos 2 ib — 2 cos ij, — 1) 






U{ib) = 


i ( 15 cos 2 ib — 10 cos ib — 


1) 




92{ib) = 


|(1 — cos it,) (15 cos 2 ib + 


10 cos if, 


-1) 


h 3 (ib) = 


i (5 cos 2 ib + 2 cos it — 1) 






hi{i b ) = 


|(15 cos 2 ij, + 10 cos if, — 


1) 





de b 
dt 

dk 
dt 

dzu 



Hi 

8 °6 



- 8°f> 



C ( h q) e b sin 2 i b sin(2w 6 ) - ^C ( b o) e c [f 2 {i b ) sm{zu b - vj c ) + g 2 {i b ) sin(zu b + w c - 2fi 6 )] + 0(x 5 ), (A12) 
siri(Slt, — Q c ) cos if, sin 2i c — ffC^ebec sin(2it,) [fi(ib) s'm(zub — zu c ) — h,4,(i b ) sin(n7t, + zo c — 2Q. b )] + 0(x°), (A13) 
^ = jCi q) [h(i b ) + § sin 2 i b cos(2wi,)] - ±|cf > (g) [f 2 {i b ) cos(w b - vj c ) + g 2 (i b ) cos( CTi) + w c - 20 6 )] + 0(x 4 ), (A14) 

(A15) 



— a n 

— ~4°t> 



(9) 



cos it, + 0(x ), 



|C*< 9) e- 1 sin(n l ,-a 



~~dt 

^ = f a Ci o) e b [f 2 (i b )sm(zo b 

di c 
dt 
dm 
~~dt 
and 

HT = ~» 

where 



- tu c ) — 32(4) sin(ci7(, + zu c — 2fl b )] + 0(x 
I sin2if, cosic + 0(x 11 '' 2 ), 



13/2s 



= |C<"/i(i 6 ) + ^ 9/2 ) 



(«)« 



_1 cos(Slt, 



o \ • o- cos2i c 7 / 2 
S2 C ) sin 2«i,— — ; h u(x 1 ) 



C, 



(9) 



«i, 



-.(9) 



In addition 



/mc\ /at 
Vm*/ \a c 

/mt\ /at 
\m, / \a c 



(o) 




and C, 



(°) 



/mt\ (Ob 



\m c J 



(5) - 



Wb 



Sib 



(A16) 
(A17) 
(A18) 

(A19) 

(A20) 
(A21) 

(A22) 



Our decision to retain or ignore each particular term is guided by numerical solutions for the relaxed state of systems with 
arbitrary relative inclinations. For example, we have omitted the quadrupole term — -jfC^'e 2 sin(2it,) sin(2ojf,) in dib/dt which 
is responsible for Kozai oscillations and has a modulation frequency of 2tjf,; in the relaxed state this term is 0(x ) and 
contributes negligibly to the dynamics, whereas the quadrupole term we have included dominates the behaviour of ib with its 



11 Note that since our TZ has the dimensions of energy, the usual equations are divided by [i b and fj, c for the rates of change of the inner 
and outer orbital elements respectively. 
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modulation frequency of cj),. Moreover, while the octopole term is 0(x 5 ), we have included it because it provides a timescale 
for the slow decay of the mutual inclination (see Section 13. 5[) . 



